############################################################################
###ELECTION HEALTH ANALYSIS ANALYSIS
###KBS MARCH 2021
###
####  
#############################################################################


####################SET WORKING DIRECTORY AND LOAD LIBRARIES##################

setwd ("C:/Users/ksmith1/Dropbox/Politics and Health/Data/Master Data and Codebooks/Working Data/PLoS Replication Posted on Dataverse") #workding director form office

library(plm)
library(foreign)
library(clusterSEs)
library(lmtest)
library(rms)
library(stargazer)
library(lmtest)
library(boot)
library(ExtremeBounds)
library(psych)

###########################LOAD DATA#################################

#data is an spps working file built from master data set

mydata <- read.spss("panel_data.sav", to.data.frame=TRUE)

str(mydata)

#just changing the pid and political interest variables to numeric

mydata$newsint <- as.numeric(mydata$newsint)
mydata$dem_ind_rep <- as.numeric(mydata$dem_ind_rep)

str(mydata)

head(mydata)
summary(mydata)



#######################################POOLED MODELS###########################################
#
# These models are used to generate the estimates reported in Table 2 of the manuscript       #
# the Hausman test for random versus fixed effects is reported for each model in Table S9     #
#                                                                                             #
###############################################################################################


################################################HEALTH SCALE MODEL################################


pooled.health <- plm(formula=health6~age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                   + newsint + political_knowledge + political_participation + TrumpVote 
                   + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.health)



# Random effects model

re.health <- plm(formula=health6 ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                 + newsint + political_knowledge + political_participation + TrumpVote 
                 + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.health)


#fixed effects


fe.health <- plm(formula=health6 ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                 + newsint + political_knowledge + political_participation + TrumpVote 
                 + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.health)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.health)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.health, pooled.health)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.health, fe.health)



#######################################################################################################



################################FULL SCALE 32 ITEM MODEL#############################################




pooled.full <- plm(formula=full_32item_health ~ age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                     + newsint + political_knowledge + political_participation + TrumpVote 
                     + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.full)



# Random effects model

re.full <- plm(formula= full_32item_health ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                 + newsint + political_knowledge + political_participation + TrumpVote 
                 + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.full)


#fixed effects


fe.full <- plm(formula=full_32item_health ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                 + newsint + political_knowledge + political_participation + TrumpVote 
                 + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.full)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.full)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.full, pooled.full)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.full, fe.full)


####################################EMOTIONAL MODELS################################################





pooled.emot <- plm(formula=Emotional_8_item ~ age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                   + newsint + political_knowledge + political_participation + TrumpVote 
                   + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.emot)



# Random effects model

re.emot <- plm(formula= Emotional_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                 + newsint + political_knowledge + political_participation + TrumpVote 
                 + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.emot)


#fixed effects


fe.emot <- plm(formula=Emotional_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.emot)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.emot)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.emot, pooled.emot)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.emot, fe.emot)





###################################SOCIAL LIFESTYLE BATTERY#####################################





pooled.soc <- plm(formula=Social_Lifestyle_8_item ~ age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                   + newsint + political_knowledge + political_participation + TrumpVote 
                   + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.soc)



# Random effects model

re.soc <- plm(formula= Social_Lifestyle_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.soc)


#fixed effects


fe.soc <- plm(formula=Social_Lifestyle_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.soc)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.soc)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.soc, pooled.soc)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.soc, fe.soc)



###########################################COMPULSION MODELS##########################################




pooled.comp <- plm(formula=Compulsion_Battery ~ age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                  + newsint + political_knowledge + political_participation + TrumpVote 
                  + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.comp)



# Random effects model

re.comp <- plm(formula= Compulsion_Battery ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
              + newsint + political_knowledge + political_participation + TrumpVote 
              + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.comp)


#fixed effects


fe.comp <- plm(formula= Compulsion_Battery ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
              + newsint + political_knowledge + political_participation + TrumpVote 
              + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.comp)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.comp)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.comp, pooled.comp)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.comp, fe.comp)




##################################10 ITEM SHORT FORM MODELS##########################################




pooled.short <- plm(formula=political_malady_additive ~ age + male + black + dem_ind_rep  ++ neg_part_score + resilience_score
                   + newsint + political_knowledge + political_participation + TrumpVote 
                   + post_election_dummy, data=mydata, model="pooling", index=c("caseid_w1", "Index1"))
summary(pooled.short)



# Random effects model

re.short <- plm(formula= political_malady_additive ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.short)


#fixed effects


fe.short <- plm(formula= political_malady_additive ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy, data=mydata, model="within", index=c("caseid_w1", "Index1"))
summary(fe.short)




# LM test for random effects versus OLS -- if significant more support for random effects over OLS

plmtest(pooled.short)

# LM test for fixed effects versus OLS -- if significant more support for FE over OLS--not reported in ms
#but shows FE preferred over OLS
pFtest(fe.short, pooled.short)


# Hausman test for fixed versus random effects model-- if significant need to go with fixed effects
phtest(re.short, fe.short)



#############################TABLE OUTPUT OF RANDOM EFFECTS MODELS#################################
#
#  This just generates a combined table of all the RE models--this is reported as Table 2 in text #



stargazer(re.full, re.health, re.emot, re.comp, re.soc, re.short, digits=2, type="text", out="re_panel_models_updated.txt")




##############################INTERACTION MODELS#####################################################
#                                                                                                   #
#   These models are identical to those reported in Table 2 with interaction terms added            #
#                                                                                                   #
#####################################################################################################

re.interact.full <- plm(formula= full_32item_health ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
               + newsint + political_knowledge + political_participation + TrumpVote 
               + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
               + polpartXelect + polknowXelect +negpartXelect + ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.full)

re.interact.health <- plm(formula= health6 ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                        + newsint + political_knowledge + political_participation + TrumpVote 
                        + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
                        + polpartXelect + polknowXelect +negpartXelect +ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.health)

re.interact.emotion <- plm(formula= Emotional_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                        + newsint + political_knowledge + political_participation + TrumpVote 
                        + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
                        + polpartXelect + polknowXelect +negpartXelect + ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.emotion)

re.interact.social <- plm(formula= Social_Lifestyle_8_item ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                        + newsint + political_knowledge + political_participation + TrumpVote 
                        + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
                        + polpartXelect + polknowXelect +negpartXelect + ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.social)


re.interact.compulsion <- plm(formula= Compulsion_Battery ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                          + newsint + political_knowledge + political_participation + TrumpVote 
                          + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
                          + polpartXelect + polknowXelect +negpartXelect +ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.compulsion)


re.interact.short <- plm(formula= political_malady_additive ~ age + male + black + dem_ind_rep + neg_part_score + resilience_score
                              + newsint + political_knowledge + political_participation + TrumpVote 
                              + post_election_dummy + ageXelect + blackXelect + maleXelect + politnXelect
                              + polpartXelect + polknowXelect +negpartXelect + ideoXelect, data=mydata, model="random", index=c("caseid_w1", "Index1"))
summary(re.interact.short)


###################################################################################################

#          This creates a table combining results of all the interaction models                   #
#           The table outputted is reported as Table S8 in the supplementary materials            #
###################################################################################################


stargazer(re.interact.full, re.interact.health, re.interact.emotion, re.interact.compulsion, re.interact.social, re.interact.short, digits=2, type="text", out="re_interact_panel_models.htm")


